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Understanding the origin of the accelerated expansion of the Universe poses one of the greatest 
challenges in physics today. Lacking a compelling fundamental theory to test, observational efforts 
are targeted at a better characterization of the underlying cause. If a new form of mass-energy, 
dark energy, is driving the acceleration, the redshift evolution of the equation of state parameter 
w(z) will hold essential clues as to its origin. To best exploit data from observations it is necessary 
to develop a robust and accurate reconstruction approach, with controlled errors, for w(z). We 
introduce a new, nonparametric method for solving the associated statistical inverse problem based 
on Gaussian Process modeling and Markov chain Monte Carlo sampling. Applying this method to 
recent supernova measurements, we reconstruct the continuous history of w out to redshift z = 1.5. 

PACS numbers: 98.80.-k, 02.50.-r 
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Little more than a decade has passed after supernova 
observations first found evidence for the accelerated ex- 
pansion of the Universe [l[ . Since confirmed by different 
probes, this remarkable discovery has been hailed as the 
harbinger of a revolution in fundamental physics and cos- 
mology. Cosmic acceleration demands completely new 
physics - it challenges basic notions of quantum theory, 
general relativity, and assumptions regarding the funda- 
mental make-up of matter. Currently, the two most pop- 
ular explanations are a dark energy, usually modeled by 
a scalar field, or a modification of general relativity on 
cosmic length scales. In the absence of a compelling can- 
didate theory to explain the observations, the target of 
current and future cosmological missions is to character- 
ize the underlying cause for the accelerated expansion. 
In the case of dark energy, the aim is to constrain the 
equation of state parameter w — p/p and its possible 
evolution. A deviation from w — const, would provide 
clues pointing to the origin of the accelerated expansion. 
(Currently, observations are consistent with a cosmolog- 
ical constant, w = —1, at the 10% level @.) 

In order to extract useful information from cosmolog- 
ical data, a reliable and robust reconstruction method 
for w(z) is crucial. Here, we introduce a nonparamet- 
ric technique based on Gaussian Process (GP) modeling 
and Markov chain Monte Carlo (MCMC) sampling, and 
apply it to supernova data. GPs extend the multivari- 
ate Gaussian distribution to function spaces, with infer- 
ence taking place in the space of functions. The defining 
property of a GP is that the vector that corresponds to 
the process at any finite collection of points follows a 
multivariate Gaussian distribution. Gaussian processes 
are elements of an infinite dimensional space, and can 
be used as the basis for a nonparametric reconstruction 
method. GPs are characterized by mean and covariance 
functions, defined by a small number of hyperparame- 
ters 0. The covariance function controls aspects such 



as roughness of the candidate functions and the length 
scales on which they can change; aside from this, their 
shapes are arbitrary. Bayesian estimation simultaneously 
evaluates the GP hyperparameters (so-called to prevent 
confusion with the parameters that define a parametric 
method) together with quantities of physical interest. 

For supernovae, the reconstruction task can be sum- 
marized as follows. The data is given in the form of the 
distance modulus hb(z) defined as: 

f i B {z)=m B - M B = 51og 10 (^^\ + 25, (1) 

where m B and M B are the apparent and absolute mag- 
nitudes. The luminosity distance cIl(z) is connected to 
the Hubble expansion rate H(z), and thus to w{z), via: 
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where h(z) = H(z)/H . Note that supernovae cannot be 
used to determine Hq in the absence of an independent 
distance measurement. Thus it remains an unknown and 
can be absorbed in a redefinition of the absolute magni- 
tude: M B = M B - 51og 10 (iT ) + 25. The H used to 
obtain A4b is not the physical value, but a certain fixed 
constant assumed in the observational analysis. For the 
dataset analyzed below [2j, Hq — 65 km/s/Mpc. Since we 
will work with n B and not with tub directly, the numeri- 
cal value for does not enter in our analysis. We allow 
for a free constant, A M in Eqn. ([T]), accounting for the 
uncertainty in the absolute calibration of the data, and 
for which we choose a uniform prior between [—0.5, 0.5]. 
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We assume spatial flatness as an "inflation prior" ; strong 
constraints on spatial flatness exist from combining cos- 
mic microwave background (CMB)and baryon acoustic 
oscillation (BAO) measurements | ]| . The prior on fl m 
is also informed by the 7-year WMAP analysis [J] for a 
wCDM model combining CMB, BAO, and Hq measure- 
ments. Since our assumptions on w are less strict than 
w = const, we broaden the nominal range by a factor 
of two, leading to a Gaussian prior with mean 0.27 and 
standard deviation of 0.04. 

Reconstruction is a classic statistical inverse problem 
for the nonlinear (smoothing) operator of Eqn. ([2]) , where 
one solves for the function w{z) and for the parameter 
£7 m , given a finite set of noisy data for cIl(z). Different 
strategies may be followed to arrive at a tractable for- 
mulation, (i) Assume a parameterized form for w(z) and 
estimate the associated parameters. This approach is the 
most commonly used currently; the parametric forms ei- 
ther assume w — const, or allow for a fixed redshift vari- 
ation such as w = wo — W\z/(1 + z), where wq and wi 
are constants Q. (ii) Pick a simple local basis repre- 
sentation for w[z) (bins, wavelets), and estimate the as- 
sociated coefficients (effectively a piecewise constant de- 
scription), using Principal Component Analysis (PCA) if 
needed, to work with eigenmodes defined as linear com- 
binations of bins @, 0] • (hi) Follow a procedure similar to 
(ii) - without PCA - but actually use (filtered) numerical 
derivatives to estimate w(z) @. (iv) Use a distribution 
over (random) functions that can represent w(z) and es- 
timate the statistical properties thereof, given observed 
data. Methods (i), (ii), and (iv) can all be carried out 
using a Bayesian viewpoint and exploring posteriors by 
MCMC methods, whereas (iii) - as presented in the liter- 
ature - represents a different class of approach to the in- 
verse problem. Taking numerical derivatives is generally 
a difficult task and a corresponding error theory seems 
hard to develop. Approach (i) can encounter difficulties 
if w(z) has a nontrivial evolution. The finite parameteri- 
zation and the specific functional form assumed can bias 
results for the temporal behavior of w(z) 9]. Methods 
(ii) and (iv) apply different philosophies - (ii) applies a 
local view of the reconstruction (z bins), whereas (iv) at- 
tempts to sample the posterior continuously in z. In a 
mild sense, the choice of a piecewise continuous represen- 
tation in (ii) - of which, w = const, is just the one-bin 
limit - can force an unphysical view of w(z), since the ac- 
tual w(z) is not piecewise constant. In contrast, method 
(iv) while fully nonpar ametric, is potentially more gen- 
eral and flexible compared to the other methods. 

Our new GP modeling-based approach is a realization 
of method (iv). It enables the identification of nontriv- 
ial redshift dependences in w(z) reliably, if they exist 
(Ref. |l0| shows examples based on simulated data) . The 
central idea is to assign prior probabilities to classes of 
functions via GPs and to take advantage of the particular 
integral structure of Eqn. ([2"]l. again using GPs. Employ- 
ing a Bayesian approach to explore posterior distribu- 
tions over the functions via MCMC we not only obtain 
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FIG. 1. Priors (red lines) and posteriors (black lines) for the 
GP hyperparameters p and k 2 . The lower left panel shows the 
distribution of the GP mean. The lower right panel shows the 
results for A M . The posteriors for different a choices are very 
similar, and we show only the results for a ~ 2 here. 

a family of continuous realizations for w(z) but at the 
same time optimize the GP model hyperparameters, in- 
formed by the actual data, comprehensively propagating 
all estimation uncertainties. 

One may wonder whether a general nonparametric re- 
construction must involve taking a second derivative of 
the data in some way. This is true only in a formal 
sense - the approach described here does not involve any 
derivatives. Instead, we invert an integral equation, ill- 
posed because the operator to be treated is a complicated 
smoothing operator involving two integrals. To make the 
problem well-behaved we make mild continuity assump- 
tions about w(z) which are justified if the origin of dark 
energy is to be described by a reasonable physical model. 

As previously noted, a GP model assumes that 
w(zi), w(z n ), for any set of redshifts z±,...,z n , follow 
a multivariate Gaussian distribution specified by mean 
and covariance functions 3] . Here we use a mean of neg- 
ative one as the prior and exponential family covariance 
functions written as (p being a numerical constant): 

K(z,z') = K 2 p |z_2 ' r - 

The hyperparameters p S (0, 1) and k, and the pa- 
rameters defining the likelihood, are determined by the 
data. The value of a <E (0, 2] influences the smooth- 
ness of the GP realizations: for a = 2, the realiza- 
tions are smooth with infinitely many derivatives (this 
covariance function corresponds to using an infinite num- 
ber of Gaussian basis functions), while a — 1 leads to 
rougher realizations suited to modeling continuous non- 
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(mean-squared)-differentiable functions. We use both 
a values in our analysis, the results being very simi- 
lar. The form of the covariance function is unrelated to 
the shape of the GP sampling functions, as determined 
by the data, so no particular behavior is assumed for 
w(z). There is no loss of generality in fixing the (sta- 
tistical ensemble) GP mean, as any variation imposed 
by the data appears in the covariance function. Fix- 
ing the mean has the advantage of improving the sta- 
bility of the MCMC (we explored other means and found 
very similar results). We stress that even though the 
mean is fixed, each GP realization will actually have a 
different mean with a spread controlled by k as shown 
in Figure [1] The constant p has a prior of Beta(6, 1) 
and k 2 has a vague Inverse Gamma prior IG(6, 2). The 
probability distribution of the Beta prior is given by 
f(x;a,/3) = T{a + p)x a - 1 {\ ~ x)^ l /[T{a)T(j3)} and for 



the IG prior by f(x; a, /3) = fi c 



T(a)- 1 exp(-/3/a;), 



with x > 0. We show the priors for p and k 2 and their 
posterior distribution in Figure [TJ 

Following the notation of Eqn. © we set up the fol- 
lowing GP for w: 

w(u) ~ GP(-l,K(u,u')). 

Recall that we have to integrate over w(u) (Eqn. [J): 

w(u) 



o 1 + w 



-du. 



(3) 



We use the fact that the integral of a GP is also a GP with 
mean and covariance dependent on the original GP 
Using that result, we set up a second GP for y(s): 

p W-^'\ a dudu' \ 
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The mean value for this GP is simply obtained by solv- 
ing the integral in Eqn. ^ for the mean value of the GP 
for w(u), here taken to be negative one. As mentioned 
earlier, even though the ensemble mean is fixed to a con- 
stant value at any z, each GP realization is not mean-zero 
(over z). We show the distribution of the mean for the 
different realizations in Figure Q] demonstrating a wide 
spread between -2 and -0.6. In addition, we used simu- 
lated data with a time- varying equation of state to check 
that the choice of the mean does not bias the results. 
We can now construct a joint GP for y(s) and w(u): 
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FIG. 2. Nonparametric reconstruction of w(z) based on GP 
modeling combined with MCMC. The upper panel uses a 
Gaussian covariance function, the bottom panel, an exponen- 
tial covariance function. Both results are very close and in 
agreement with a cosmological constant (black dashed line). 
The dark blue shaded region indicates the 68% confidence 
level, while the light blue region extends it to 95%. 



The mean for y(s) given w(u) can be found through the 
following relation: 

(y(s)\w(u)) = - ln(I + a) + S 12 E 22 1 [w(u) - (-1)] . 

Now only the outer integral is left to be solved for in 
Eqn. ([2]), and this can be computed by standard numer- 
ical methods. (Note that the computationally expensive 
double integral for En as defined in Eqn. ((4]) does not 
need to be performed.) The GP prior can now be com- 
bined with a likelihood function to obtain a posterior 
that can be sampled by MCMC methods. Details of our 
GP-based MCMC implementation are provided in the 
supplementary material j7l| . 

For our specific analysis we focus on a recent compos- 
ite supernova dataset, provided by Hicken et al. B|. This 
dataset combines the so-called Union dataset [12| with 
new measurements of low redshift supernovae to form 
the Constitution set. The dataset has been analyzed in 
Ref. 0] using different light curve fitters for the supernova 
light curves; our analysis uses results from the SALT fit- 
ter (Table I in Ref. [2], which includes estimates for the 
error for the distance modulus /is - the tables contain 
what is referred to in Ref. @ as the "minimal cut" ) . 
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Our final results for w(z) are shown in Figure [2j The 
upper panel shows the results from a GP model with a 
Gaussian covariance function (a ~ 2) while the results 
in the lower panel are based on an exponential covari- 
ance function (a = 1). The results are very similar, 
the Gaussian covariance function leading to a slightly 
smoother prediction. The mean value of w(z) is very 
close to -1 at redshifts close to zero and rises slightly 
at redshift z — 1.5. Within our estimated errors, the re- 
sults are consistent with a cosmological constant w = — 1. 
Note, however, that realizations of w(z) with nontriv- 
ial z-dependence are not excluded; as observations im- 
prove the allowed range of variability will be further con- 
strained. In Rcf. ]2| a combined analysis of supernova 
data and baryon acoustic oscillation measurements is car- 
ried out. Assuming w — const., the SALT-based dataset 
yields w = — 0.987iQ'og§ consistent with our findings. 

To summarize, we have presented a new, nonparamet- 
ric reconstruction technique for the dark energy equation 
of state and applied it to current supernova observations. 
The GP-based method can be used to determine the most 
probable behavior of w(z) and to infer how likely a tar- 
get trajectory is given the current data. Thus it can 
be used to accept or reject classes of w(z) models. The 
method allows adjusting of smoothness assumptions re- 
garding w(z); priors on the GP hyperparameters control 
the allowed arbitrariness (e.g., degree of differentiabil- 
ity). Robustness of the results obtained can be checked 



by varying these priors. Our results for w(z) are consis- 
tent with a cosmological constant, with no evidence for a 
systematic mean evolution in w with redshift, although 
variations within our error limits cannot be ruled out. We 
have carried out careful tests to ensure that our choices of 
priors and hyperparameters do not alter the results. Our 
method possesses several advantages: it avoids artificial 
biases due to restricted parametric assumptions for w(z), 
it does not lose information about the data by smooth- 
ing it, and it does not introduce arbitrariness (and lack of 
error control) in reconstruction by representing the data 
using a certain number of bins, or cutting off information 
by using a restricted set of basis functions to represent 
the data. The technique can be easily extended to fold 
in data from CMB and BAO observations; work in this 
direction is currently in progress. The GP-based MCMC 
procedure can be integrated within supernova analysis 
frameworks, e.g., SNANA [l3{ as a cosmology fitter, fol- 
lowing the general methodology presented in Ref. [14j . 
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